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. ■ Abstract. - 

A model of a finite cylindrical ion channel through a phospholipid membrane of width L separating 
two electrolyte reservoirs is studied. Analytical solution of the Poisson equation is obtained for an ar- 
bitrary distribution of ions inside the trans-membrane pore. The solution is asymptotically exact in the 
limit of large ionic strength of electrolyte on the two sides of membrane. However, even for physiological 
concentrations of electrolyte, the electrostatic barrier sizes found using the theory are in excellent agree- 
ment with the numerical solution of the Poisson equation. The analytical solution is used to calculate 
O , the electrostatic potential energy profiles for pores containing charged protein residues. Availability of a 

semi-exact interionic potential should greatly facilitate the study of ionic transport through nanopores and 
^vq ion channels. 
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Ion channels are water filled holes which facilitate exchange of electrolyte between the exterior and 
interior of a cell. Pores are formed by specific proteins embedded into the phospholipid membrane [ 1 ] . 
Depending on the conformation of the protein, the pore can be open or closed. When open, the protein 
is very specific to the kind of ions that it allows to pass through the channel [2, 3]. In order to function 
properly the channel has to conduct thousands of ions in a period of few milliseconds. Considering 
that the channel passes through a phospholipid membrane which has a very low dielectric constant 
and is very narrow, producing high energetic penalties for ions entering the nonopores, it is fascinating 
to contemplate how Natures manages to perform this amazing task. In fact, as long ago as 1969, 
Parsegian observed that for an infinitely long cylindrical channel [4] of radius a = 3 A, the electrostatic 
barrier is over 16fcsT, which should completely suppress any ionic flow [5], Later numerical work 
by Levitt [6], Jordan [7] and others demonstrated that for more realistic finite channels the barrier is 
dramatically reduced. For example, for a channel of length L — 25 A and radius a = 3 A, the barrier 
5_i ■ is about 6fcsT, which although still quite large, should allow ionic conductivity. Recently the study 
of ion channels has expanded to other parts of applied physics. Water filled nanopores are introduced 
into silicon oxide films, polymer membranes, etc [8, 9]. In all of these cases the dielectric constant of 
the interior of a nanopore greatly exceeds that of the surrounding media. 

To quantitatively study the conductance of a nanopore one has three options: the all atom molecular 
dynamics simulation (MD) [3]; the Brownian dynamics simulation (BD) [10, 11] with implicit water 
treated as a uniform dielectric continuum; or the mean-field Poisson-Nernst-Planck theory (PNP) [12] 
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which treats both water and ions implicitly. While clearly the most accurate, MD simulations are 
computationally very expensive [13]. Brownian dynamics is significantly faster than MD, but because 
of the dielectric discontinuities across the various interfaces a new solution of the Poisson equation is 
required for each new configuration of ions inside the pore. The simplest approach to study the ionic 
conduction is based on the PNP theory [12], This combines the continuity equation with the Poisson 
equation and Ohm's and Fick's laws. PNP is intrinsically mean-field and is, therefore, bound to fail 
when ionic correlations become important. This has been well studied for its static version — the 
Poisson-Boltzmann equation, which is known to break down for aqueous electrolytes with multivalent 
ions and also for monovalent electrolytes in low dielectric solvents [14, 15]. For narrow channels, 
the cylindrical geometry, combined with the field confinement, results in a pseudo one dimensional 
potential of very long range [16, 17]. Under these conditions the correlational effects dominate, and the 
mean-field approximation fails [14]. Indeed recent comparison between the BD and the PNP showed 
that PNP breaks down when the pore radius is smaller than about two Debye lengths [10, 1 1]. At the 
moment, therefore, it appears that a semi-continuum (implicit solvent) Brownian dynamics simulation 
is the best compromise between the cost and accuracy [13,18,19] for narrow pores. Unfortunately even 
this, simplified strategy demands a tremendous computational effort. The difficulty is that BD requires 
a new solution of the Poisson partial differential equation at each time step. This can be partially 
overcome by using lookup tables [11] and variational methods [20], but still requires a supercomputer. 
If the interaction potential between the ions inside the channel would be known, the simulation could 
proceed orders of magnitude faster. However, up to now the only exact solution to the Poisson equation 
in a cylindrical geometry was for the case of an infinitely long pore [4, 5, 16]. In this letter we shall 
provide another exact solution, but now for a finite trans-membrane channel. 

We shall work in the context of a primitive model of electrolyte and membrane. The membrane 
will be modeled as a uniform dielectric slab of width L located between z — and z — L. The 
dielectric constant of the membrane and the channel forming protein is taken to be e p w 2. On both 
sides of the membrane there is an electrolyte solution composed of point-like ions and characterized 
by the inverse Debye length n. A channel is a cylindrical hole of radius a and length L filled with 
water. As is usual for continuum electrostatic models [13], we shall take the dielectric constant of 
water inside and outside the channel to be the same, e w ss 80. It is convenient to set up a cylindrical 
coordinate system (z, p, <j>) with the origin located at the center of the channel at z = 0. Suppose that 
an ion is located at an arbitrary position x' inside the channel. The electrostatic potential ip(z, p, <p; x') 
inside the channel and membrane satisfies the Laplace equation 

VV = -^(x-x'). (1) 

For z > L and z < 0, <^(x;x') satisfies the linearized Poisson-Boltzmann or the Debye-Hiickel 
equation [14] 

VV = K 2 ip . (2) 

The inverse Debye length is related to the ionic strength / of electrolyte, ^ = n = \/8tt\bI, where 
\b = q 2 /eksT is the Bjerrum length and / = (a 2 c a + ac a + 2c)/2. Here c a is the concentration 
of a : 1 valent electrolyte and c is the concentration of 1 : 1 electrolyte. All the usual boundary 
conditions must be enforced: the potential must vanish at infinity and be continuous across all the 
interfaces; the tangential component of the electric field and the normal component of the electric 
displacement must be continuous across all the interfaces. These boundary conditions guaranty the 
uniqueness of solution. Unfortunately, even this relatively simple geometry can not, in general, be 
solved exactly. We observe, however, that an exact solution is possible in the limit that k — ► oo. In 
this special case the system of differential equations becomes separable. Our, strategy then will be to 
solve exactly this asymptomatic problem and then extend the solution to finite values of Debye length. 
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We start by making the following fundamental observation. The condition n — > oo signifies that 
electrolyte perfectly screens any electric field — the Debye length is zero. This, combined with the 
boundary condition — electrostatic potential must vanishes at infinity — implies that in this limit 
<p(0, p, cj>; x') = ip(L, p, (f>; x') = 0, for any position x' of an ion inside the pore. This is a dramatic 
simplification. Now it is no longer necessary to solve the Debye-Hiickel equation, but only the Poisson 
equation with a perfect grounded conductor boundary conditions at z = and z = L. To proceed we 
expand the 8{z — z 1 ) in eigenfunctions of the differential operator 



dz 2 



+ K^Pn = , (3) 



satisfying the perfect conductor boundary condition. The normalized eigenfunctions are ip n (z) = 
y/2/Lsm{k n z), with k n — nir/L. The Sturm-Liouville nature of the differential Eq. guaranties 
us that 



Similarly, 



Next we write 



S(z - z') = — 2_^ sm(k n z) sin(fc„z') . (4) 



_. oo 

— J2 e """" rt • (5) 

m= — oo 



^(x,x') = ^ T £ JT e im ^-^8m(k n z)sia(k n z')g nm (p,p'). (6) 

n—1 m— — oo 

Substituting this into Eq.Q we find that the Green function g nm {p, p') satisfies the modified Bessel 
equation 

1 d dg nm 2 to 2 4tt 



solution of which can be found using the usual techniques [4, 21]. We obtain, 

9mn(p,p') = ^Im(knP<)[K m (knP>) + lmnI m {k n P>)] (8) 

where p> and p < are the larger and the smaller of the set (p, p') and 

K m (k n a)K' m (k n a)(e p - e w ) 

^w i- m\kn a ) "■ m\knOi) €pl rn yk n a,)l\ rn yk n o 1 ) 

Here I m , K m , I' m , K' m are the modified Bessel functions of the first and second kind and their deriva- 
tives, respectively. Eqn.(|8j is valid for p > < a. When p > a, 

_ 4ttc w K m (k n p)I m (k n p') 

k n a E w l rn yk n a)l\ rn yK n a) £pl m (k n a)K nl (k n a) 

Eans. d6l8t . and dlOt are exact for an ion inside a pore with perfect conductor boundary conditions at 
z = and z = L. If the ion is located on the axis of symmetry, z' = zq, p' = 0, only m = term in 
Eqn.(|6j survives, and the electrostatic potential inside the channel at position z, p takes a particularly 
simple form, <pi n (z, p; z ) = ipi(z, p; z ) + ip 2 (z, p; z ), where 

<Pi(z,P'>Zo) = — fy^sm(k n z)sm(k n z )K (k n p) , (11) 

e w L ^ 

n=l 
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and 

Aq(e w - e p ) ^ K Q (k n a)Ki(k n a)Io(k n p)sm(k n z)sm(k n zo) 

(p 2 {z,p;z Q ) = > — — . (12) 

e™^ ^ e w li(k n a)K {k n a) + eplo(k n a)Ki{k n a) 

Eqs. i ll 1112b are exact in the « — * oo limit. To see how these equations can be extended to finite values 
of k, it is important to first understand their physical meaning. Potential if 2 is mostly the result of the 
charge induced on the interface between the high dielectric aqueous interior of the pore and the low 
dielectric membrane. We expect that this term will be affected very little by the precise value of the 
Debye length of the surrounding electrolyte solution. 

The potential ipx contains the contribution from the ion located at zo and from the induced charge 
on the pore/electrolyte and the membrane/electrolyte interfaces. It will, therefore, strongly depend 
on the precise value of n. Furthermore, we observe that Eq. dl It is exactly the potential produced 
by a charge q located inside an infinite slab of water of width L bounded by two grounded perfectly 
conducting planes [21]. This key observation allows us to explicitly resum the series in Eq. i ll It . 
However, it is possible to do even better, and now enforce the exact boundary condition, namely that 
for z < and z > L the electrostatic potential must satisfy the Debye-Huckel equation (|2j. Using the 
Bessel J representation of the delta function one can constructs the Green function [22] which satisfies 
all the boundary conditions for the slab geometry and has the required symmetry property [2 1 ] between 
the source and the observation points. We then find 

f°° Ja(kp){a 2 (k)e klz - z ^- 2kL +a(k)P(k)[e' k{z+Zo) +e k{z+z °^ 2kL }+P 2 (k)e- klz - Zo1 } 

viO,p; 2 o)= / dfc 



/3 2 (fc)-a 2 (fc)exp(-2fcL) 

(13) 
where a(k) = [k - \/k 2 + n 2 ]/2k, f3(k) = [k + \Jk 2 + n 2 ]/2k, and J (x) is the Bessel function of 
first kind and order zero. Eq. Jl 31 provides an analytic continuation of Eq. Jill into finite k parameter 
space. It can be checked explicitly that in the limit k — > 00, Eq. Jl 3i exactly sums the series in 
Eq. ( Ill> . Finally, for the region p > a the electrostatic potential is 

4q ^ 1 K Q (k n p)sm(k n z)sm(k n z ) 

(p out {z,p;z ) = — > — — — — — - — — — — — — — -. (14) 

=1 ™ a € wh{Kna)Ko(kna) + e p 7o(K„a)Ai(fc„a) 

Eq J14l > is exact only for the perfect conductor boundary conditions, however, the huge jump in the di- 
electric constant going from the membrane's interior to the aqueous electrolyte will leave (f ou t mostly 
unaffected even for finite values of n. 

If the channel contains N ions and charged protein residues their interaction energy is given by 



1 N 



where qi is the charge of ion/residue % and tp 3 is the electrostatic potential produced by the ion/residue 
j at the position of ion/residue i. Similarly the electrostatic barrier that an ion feels as it moves through 
a charge free channel is [14], 



U(z) = - lim 
1 p^O 



tp(z,p;z) 



Q 






The last term in Eq. d!6t is the electrostatic "solvation" energy that a point-like ion looses as it 
moves from the bulk electrolyte into the interior or a pore. This energy can be calculated using the 
Debye-Hiickel theory and is equivalent to the excess chemical potential resulting from the screening 
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of ionic electric field by the surrounding electrolyte [14]. The limit in Eq. dl6l is easily obtained by 
noting that 

- = / J (kp)dk . (17) 

P Jo 

We are now in a position to explore some of the quantitative consequences of the current theory. 
In Fig. we first plot the potential energy barrier for an ion of charge q moving through a channel 
of L = 35 A and a = 3 A and various external electrolyte concentrations. For the physiological 




Fie. 1 



Fig. 1 - Electrostatic potential barriers for an ion of charge q moving along the axis of symmetry through a channel 
of L = 35 A and a — 3 A . External electrolyte concentration from the bottom up is 0.15, 1, 2 and 3M. 

Fig. 2 - Electrostatic potential barrier for an ion of charge q moving through a pore of L = 25 A and a — 2.5 
A . External electrolyte concentration is 2.5M. Symbols are the result of a numerical integration of the system of 
Poisson-Boltzmann (bulk electrolyte) and Poisson equation (pore interior) from ref. [23] 



salt concentration (150mM) we find the barrier height to be 8.13fcsT. Using numerical solution of the 
Poisson equation, Levitt obtained a barrier of 8.48/csT. Some of the difference between the two values 
can be attributed to the fact that in numerical calculations presence of external electrolyte was not taken 
into account. As the length of the channel increases, the role of external electrolyte becomes relatively 
less important. Indeed for a channel of L = 50 A and a = 2 A, we obtain a barrier of 18.65fcsT, 
while the Levitt's numerical solution produced 17.2ksT [6] and Jordan's 18.6fcsT [7]. The only 
numerical work known to us that explicitly takes into presence existence of external electrolyte is 
ref. [23]. The authors of that paper numerically solved the non-linear Poisson-Boltzmann equation for 
external electrolyte and the Poisson equation for the interior of the channel. For a pore of L = 25 
and a — 2.5 A, and electrolyte concentration of 2.5M, they find a barrier of 9.5fcsT, while we obtain 
9.8/csT. In Fig. (|2j we compare the full electrostatic energy barrier obtained from the numerical 
solution with our analytical results. The agreement, once again, is quite good. 

Availability of a semi-exact interaction potential allows us to easily explore the potential energy 
landscape $ = V + U of an ion of charge q moving through a channel which also contains some 
fixed charged protein residues. For example, consider a channel of L = 35 A and a = 3 A and 
suppose that there is one protein residues of charge — q, embedded into the surface of the channel at 
(z = L/2, p — a, (j>). In Fig.[3]we show the potential energy profile for an ion moving along the axis 
of symmetry through such a channel. Instead of a potential barrier, this ion encounters a potential well 
of depth more than lO/c^T! It will, therefore, find it extremely difficult to pass through such a channel. 
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Fig. 3 - Potential energy profile for an ion of charge q moving along the axis of symmetry through a channel of 
L = 35 A and a = 3 A containing a protein residue of charge — q embedded into its wall (p = a) at z = L/2 
(solid line); containing two charged residues at z — 5 and z=30 A (dotted line); containing three residues at 
z — 5, 17.5 and z — 30 A (dashed line). External electrolyte concentration is 150mM. 

Fig. 4 - Potential energy profile for an ion of charge q moving along the axis of symmetry through a pore of 
L = 35 A and a = 3 A, containing two charged residues hidden in the membrane's interior (p = 6 A) located at 
z — 5 and z = 30 A. External electrolyte concentration 150mM. 



Now, suppose that two charged residues are embedded into the channel wall, one close to the 
entrance of the channel at z = 5 A and another close to its exit at z = 30 A, both at p — 3 A. 
The electrostatic potential, now develops a double-well structure, Fig. [3] Each minimum is relatively 
less deep than in a channel with only one central residue. One might then suppose that adding more 
charged residues will diminish the depth of the wells even further. This, however, is not the case. 
In Fig|5]we also show the potential energy landscape of a channel containing three uniformly spaced 
residues. Evidently instead of decreasing the depth of the potential well, it has dramatically increased! 
There is, however, a mechanism which Nature can use to diminish the depth of potential wells — 
hide the charged residues in the membrane's (or protein's) hydrophobic interior [24]. In Fig. |4] we 
plot the potential energy profile for the same channel as in Fig. [3] but with the two charged residues 
hidden in the membrane's hydrophobic interior at p = 6 A. In this case the deep potential well is 
replaced by a shallow binding site followed by an activation barrier of only 2.5fceT. This can be 
easily overcome by an external electric field or a chemical potential gradient. We find that there is 
an optimum location for hiding charged residues in order to produce the smallest barrier for channel 
penetration. The formalism developed above allows us to easily explore all the parameter space in 
order to find this optimum position. 

To conclude, we have presented an analytically solvable model of electrostatics inside an ion chan- 
nel. The solution found is exact in the limit of large electrolyte concentrations. However, comparison 
with the numerical work shows that it remains valid even at intermediate and low electrolyte con- 
centrations. The analytical solution can be used to dramatically speed up the Brownian dynamics 
simulations of ionic transport through cylindrical pores. The biological and structural information can 
be partially taken into account through a proper placement of charged protein residues. Furthermore, 
even if a more detailed atomistic molecular dynamics simulation is necessary, availability of a rapid 
Brownian dynamics model can serve for a initial exploration of the parameter space. 

In this work water inside the pore has been treated as a uniform dielectric continuum identical to 



Yan Levin: ELECTROSTATICS OF IONS INSIDE THE NANOPORES AND TRANS-MEMBRANE CHANNELS7 



the bulk. While this is acceptably for wider pores, the approximation will certainly fail for very narrow 
pores such as gramicidin. To properly account for the polarization of water in this geometry one must 
go beyond the continuum dielectric approximation [25]. Until now, the only option for these cases 
was to perform all atom molecular dynamics simulations. The current work suggests that another way 
might be possible. The continuum description with e w = 80 and e p — 2 can be used for the bulk 
water outside the channel, for the membrane, and for the trans-membrane protein, while inside the 
channel (now with ej„ = 1) one might try to obtain an accurate analytical electrostatic potential. The 
Coulomb interactions between the water molecules inside the channel could then be treated explicitly, 
without any need for continuum dielectric approximation. This would then allow to perform very fast 
molecular dynamics simulations of ionic transport, free of the drawbacks associated with the implicit 
solvent models. The work in this direction is now in progress. 

This work was supported in part by Conselho Nacional de Desenvolvimento Cientifico e Tec- 
nologico (CNPq). 
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